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ABSTRACT 

We present a new algorithm for identifying the substructure within simulated dark 
matter haloes. The method is an extension of that proposed by Tormen et al. (2004) 
and Giocoli et al. (2008a), which identifies a subhalo as a group of self-bound particles 
that prior to being accreted by the main progenitor of the host halo belonged to one 
and the same progenitor halo (hereafter 'satellite'). However, this definition does not 
account for the fact that these satellite haloes themselves may also have substructure, 
which thus gives rise to sub-subhaloes, etc. Our new algorithm identifies substruc- 
tures at all levels of this hierarchy, and we use it to determine the mass function of 
all substructure (counting sub-haloes, sub-subhaloes, etc.). On average, haloes which 
formed more recently tend to have a larger mass fraction in substructure and to be 
less concentrated than average haloes of the same mass. We provide quantitative fits 
to these correlations. Even though our algorithm is very different from that of Gao 
et al. (2004), we too find that the subhalo mass function per unit mass at redshift 
z = is universal. This universality extends to any redshift only if one accounts for 
the fact that host haloes of a given mass are less concentrated at higher redshifts, 
and concentration and substructure abundance are anti-correlated. This universality 
allows a simple parametrization of the subhalo mass function integrated over all host 
halo masses, at any given time. We provide analytic fits to this function which should 
be useful in halo model analyses which equate galaxies with halo substructure when 
interpreting clustering in large sky surveys. Finally, we discuss systematic differences 
in the subhalo mass function that arise from different definitions of (host) halo mass. 

Key words: galaxies: halo - cosmology: theory - dark matter - methods: numerical 
simulations - galaxies: interactions 



1 INTRODUCTION 

In the standard scenario of structure formation, galaxies are 
surrounded by extended dark matter haloes, which form by 
gravitational instability seeded by some initial density fluc- 
tuation field 8(x, z). According to the simple spherical col- 
lapse model (e.g. Peebles 1980), a halo of mass M collapses 
at a redshift z when the linear density fluctuation field - 
smoothed on a scale R ~ M 1 ^ 3 - first exceeds some criti- 
cal threshold (Press & Schechter 1974; Bond et al. 1991; 
Lacey&Cole 1993; Sheth & Tormen 2002). Small systems 
collapse at higher redshifts when the universe was denser, 
and then merge hierarchically to form larger haloes. Galaxy 
clusters are at the top of this hierarchy: they represent the 
most massive virialized structures in the present-day uni- 
verse, and may host thousands of galaxies. Recent studies, 
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using high resolution numerical simulations (e.g. Moore et 
al. 1998; Tormen et al. 1998; Springel et al. 2001), have 
shown that the cores of subhaloes accreted along the host 
merging history may survive until the present-time to form 
the so-called substructure population. 

Different studies of galaxy formation and evolution 
have attempted to correlate such substructures with satel- 
lite galaxies, with conflicting results. TV-Body simulations 
of galaxy-size haloes seem to predict more substructures 
than observed (Moore et al. 1999; Stoehr et al. 2003). A 
number of different solutions to this 'substructure problem' 
have been suggested, among others early cosmic reioniza- 
tion (Susa & Umemura 2004), mass loss and gas stripping 
(Kravtsov et al. 2004; Maccio et al. 2009), and a down- 
turn in the power spectrum at small scales (Kamionkowski 
& Liddle 2000). 

From the observational point of view, recent kinematic 
analyses of Milky-Way satellites (Simon & Geha 2007; 
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Walker et al. 2009) have shown that, if galaxy formation 
in low-mass dark matter haloes is strongly suppressed af- 
ter reionization, the circular velocity function of simulated 
CDM subhaloes can be brought into approximate agreement 
with the observed circular velocity function of the Milky 
Way satellite galaxies (e.g., Koposov et al. 2009). Ishiyama 
et al. (2008, 2009) claim instead that there is no missing 
satellites problem if one considers the scatter in the sub- 
structure mass function due to the different formation his- 
tory of the haloes. The idea here is that the Milky- Way and 
Andromeda haloes reside within an underdense region, so 
they accreted their substructures at above average redshifts 
(compared to other objects of the same mass today), and 
so they are expected to possess fewer satellites (see also van 
den Bosch et al. 2005). 

Knowledge of the substructure population is important, 
not only for galaxy formation studies, but also for estimating 
the expected observable 7-ray flux from dark matter parti- 
cle annihilations within the Milky Way. The expected flux 
depends strongly on the substructure population, its spatial 
distribution within the halo, and on the DM particle cross 
section (Pieri et al. 2008; Diemand et al. 2007; Giocoli et 
al. 2008b). 

The observed thickness of stellar disks in spiral galax- 
ies is another imprint of the substructure population of the 
host haloes. Semi-analytical models (Benson et al. 2004) and 
numerical simulations (Kazantzidis et al. 2008, 2009) show 
that merging events in the central region of the halo are re- 
sponsible for disk thickening. Hence, the substructure mass 
function, its redshift evolution and the satellite accretion 
rate all represent key ingredients for a more complete model 
of the disk structure of spiral galaxies, not to mention flares, 
bars, and faint filamentary structures above the disk plane. 

Different algorithms have been proposed to identify sub- 
structures in simulated dark matter haloes. All these algo- 
rithms yield subhalo mass functions that resemble a power- 
law at the low mass end with a logarithmic slope in the 
range from -0.8 to -1 (Gill et al. 2004a; Gao et al. 2004; 
De Lucia et al. 2004; Shaw et al. 2006; Giocoli et al. 2008a; 
Madau et al. 2008; Wetzel et al. 2009a), both at redshift 
z — and at higher redshifts. In this work we highlight some 
differences among four such methods and present a new sub- 
halo finder which we believe is useful for studies of galaxy 
formation and evolution in that it is directly linked to the 
merger tree of the host halo. 

In this paper we focus on the study of the "dark" sub- 
halo population and on its redshift evolution, to help clarify 
the role of dark matter and its dynamics in the structure for- 
mation history. In Section 2 we describe several substructure 
identification algorithms already in the literature, including 
the one we will be modifying. Section 3 describes the cos- 
mological iV-Body simulation we use, and the required post- 
processing. Our new algorithm to identify substructures, and 
an analytic fit to the resulting subhalo mass function is pre- 
sented in Section I. Section 5 discusses how the scatter in 
the substructure mass fraction at fixed mass and at differ- 
ent redshifts correlates with the formation history, and other 
structural parameters of the halo. A summary and conclu- 
sions are given in Section 7. In Appendix A we compare the 
results of this paper with two previous ones: Giocoli et al. 

(2008a) and Gao et al. (2004) that studied the substruc- 
tures population on the same cosmological simulation but 



with different algorithms. In Appendix B we propose some 
fitting functions useful to rescale the virial radius and mass 
when different definitions of the enclosed virial density are 
adopted. 



2 BACKGROUND 

To date, several different algorithms have been proposed to 
identify substructures within simulated dark matter haloes, 
and different definitions have been adopted by several au- 
thors. 

• Ghigna et al. (2000) used an algorithm called SKID. 
At each simulation snapshot, SKID estimates the density of 
each particle using a cubic spline kernel; each particle is then 
moved along the density gradients until it oscillates around 
a local density maximum. Particles are then linked using a 
friends-of-friends (FoF) algorithm (Davis et al. 1985), and 
the groups thus found are pruned iteratively to retain only 
self-bound particles. However, to obtain the full substruc- 
ture mass SKID requires a user-defined choice for the link- 
ing length to be used by the FoF algorithm. Also Weller et 
al. (2005) implemented an analogous algorithm to identify 
substructures by moving particles up density gradients and 
identifying groups as all particles reaching the same density 
maximum. In addition, they also implemented a method to 
identify a hierarchy of subhaloes within subhaloes. Shaw et 
al. (2007) improved the algorithm of determining the ener- 
getically bound components of a subhalo taking into account 
also all the forces, both internal and external, exerted on the 
subhalo. 

• Springel et al. (2001) developed SUBFIND, which de- 
fines 'subhaloes' as locally overdense, self-bound particle 
groups. SUBFIND runs on individual simulation snapshots, 
but can afterwards reconstruct the full merger tree of each 
subclump, by using the subhalo information from previous 
snapshots (Springel et al. 2005; Croton et al. 2006; Harker 
et al. 2006). The density of each particle (assumed as tracers 
of the three-dimensional dark matter density field) is esti- 
mated using an SPH-fashion scheme: the local smoothing 
scale is set to the distance of the Ndens nearest neighbor, 
and the density estimated by kernel interpolation over these 
neighbors. Locally overdense regions are identified by im- 
itating such a lowering of a global density threshold. Any 
locally overdense region enclosed by an isodensity contour 
that traverses a saddle point is considered as a substruc- 
ture candidate. All subhalo candidates are also examined 
and unbound particles are removed and redistributed. This 
algorithm has been used and tested by different authors Gao 
et al. (2004); De Lucia et al. (2004); Springel et al. (2005). 
In particular SUBFIND has been run on the largest existing 
simulation of the Milky Way (Springel et al. 2008a; Springel 
et al. 2008b). 

• Gill et al. (2004a) and Gill et al. (2004b) used MLAPM 
(Kncbe et al. 2001), an adaptive mesh refinement code 
for cosmological simulations, to locate substructures in host 
haloes. To identify haloes and subhaloes at each simulation 
output, they first build a list of potential centers by stor- 
ing the centroid of the densest grid point at the end of each 
grid tree's "branch" . Assuming that each density peak in the 
adaptive grid of MLAPM corresponds to the centre of a halo 
or subhalo, they step outwards in (logarithmically spaced) 
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radial bins until the mean enclosed overdensity reaches the 
critical virial value A v i r (z) (Eke et al. 1996; Bryan & Nor- 
man 1998), thus defining i? v ir- However, subhaloes do not 
extend out to their original virial radius, but rather to some 
truncation radius -Rtrunc, where an upturn in the radial den- 
sity profile is detected. This rise is encountered because sub- 
structures are embedded in the background of the host halo. 

• Tormen et al. (2004) and Giocoli et al. (2008a) devel- 
oped and optimized the code SURV, which identifies sub- 
haloes within the virial radius of the final host by follow- 
ing their progenitors from the time they were first accreted 
by the host main progenitor. Hence, SURV differs from the 
methods discussed above in that it uses prior information 
based on the merger tree of a host halo at redshift zo to 
identify its subhaloes. For each progenitor halo p that has 
been accreted onto the main progenitor, SURV identifies a 
subhalo as the subset of particles that belonged to p at the 
moment of accretion (z acc > zq) and are still part of a self- 
bound entity at zq within the corresponding tidal radius. 
For each subhalo thus identified, SURV stores the accretion 
time, Zacc, the original virial mass (i.e., the mass at z acc ), as 
well as the evolution of its orbital parameters after accretion. 
In this paper we extend SURV by following all branches of 
the merger tree (rather than just the main branch), in or- 
der to reconstruct the full hierarchy of substructure down to 
the mass resolution of the simulation (i.e., we aim to identify 
sub-subhaloes, sub-sub-subhaloes, etc). 



3 THE SIMULATION 

We used the data from GIF2 (Gao et al. 2004), a 
cosmological TV-Body simulation which is available at 
http://www.mpa-garching.mpg.de/Virgo. The simulation 
followed the evolution of 400 3 dark matter particles, each 
of mass 1.73 x 10 9 h _1 M , in a periodic cube of comov- 
ing side 110/i _1 Mpc within which the background cosmol- 
ogy was ACDM with (O m , cr 8 , h, Q. b h 2 ) = (0.3, 0.9, 0.7, 
0.0196). We made use of 50 snapshots, mostly logarithmi- 
cally spaced in time between z = 10 and z = 0. These are 
sufficiently closely spaced in time that one can reconstruct 
accurate halo and subhalo merger trees (Tormen et al. 2004; 
Giocoli et al. 2008a). See Gao et al. (2004) and Giocoli et 
al. (2008a) for more details about the GIF2 simulation and 
the post-processing. 



3.1 Merger Trees and Substructures 

For each simulation snapshot, haloes are identified using a 
spherical overdensity criterion: first we determine the local 
dark matter density at the position of each particle, i, by 
calculating the distance to the tenth closest neighbor, d^io- 
We assign to each particle a local density p^dm oc d~^ , sort 
particles in density and take as centre of the first halo the po- 
sition of the densest particle. We then grow a sphere around 
this centre, and stop when the mean density within the 
sphere falls below A v j r (z)p cr it(.z) as dictated by the spheri- 
cal collapse model. For a flat universe, this is approximated 
to better than one percent by 

Avir(z) = 9tt 2 (l + n m (zf -a[l- ^(2)]) (1) 




Figure 1. Schematic representation of the merger tree of a dark 
matter halo along discrete time steps. The dark-grey haloes on the 
left represents the evolution of the main halo progenitor, which 
accretes 'satellite' haloes A, B, C and D that give rise to a pop- 
ulation of subhaloes' 1 ) . System D has accreted satellite haloes 
itself (a, b and c), before it was accreted by the host. Those that 
survive (b and c) give rise to a population of subhaloes' 2 ) . 



with a = 0.7076 and f3 = 0.4403 (Stoehr 1999). For the 
GIF2 cosmology, A v i r = 97 times the critical density at 
z = 0, and it increases with redshift, asymptoting to 187r 2 at 
early times. In Appendix B we present fitting functions that 
allow halo radius, mass and concentration to be converted 
to those appropriate for other values of A v i r . 

At this point we assign all particles within the sphere 
to the newly formed halo, and remove them from the global 
list. We take the centre of the next halo at the position of 
the densest particle among the remaining ones, and grow a 
second sphere. We continue in this manner until all particles 
are screened. We include in our catalogue only haloes with 
at least 10 particles within the virial radius; particles not 
ending up in any halo are considered as 'field' or 'dust' par- 
ticles. Hereafter we will refer to the virial mass of a (host) 
halo thus identified at redshift z as M z . 

To construct merger trees we proceed as follows: for 
each halo M Zi , identified in the simulation snapshot cor- 
responding to redshift z;, we define the progenitors at the 
previous output time (at Zi+i = Zi + Az) as being all haloes 
containing at least one particle that at Zi will be in M Zi . The 
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observation redshift zq 


11.5-12 


12-12.5 


12.5-13 


13-13.5 


13.5-14 


14-14.5 


> 14.5 





8305 (10230) 


3349 (3897) 


1186 (1346) 


461 (503) 


127 (141) 


35 (36) 


4 GO 


0.5 


9347 (11252) 


3544 (39^0) 


1244 (J352) 


394 (411) 


94 ( 94) 


21 (21) 


2 (2) 


1 


9574 (11115) 


3455 (3808) 


1095 (1157) 


279 


57 (57) 


2 (2) 


1 (J) 


2 


8465 (9365) 


2461 (25S6) 


593 (605) 


98 (JOi) 


5 (5) 


2 (2) 




4 


2847 (2S7S) 


427 


35 (35) 


- (-0 









Table 1. The number of host haloes in each mass bin and at each 'observation redshift' zq. Note that we only consider host haloes whose 
mass never exceeds the mass at z = zq by more than 10%. The number in parenthesis is the total total number of haloes in each mass 
bin, including those that do not meet this criterion. 



progenitor providing the largest mass contribution to the 
parent halo is termed the main progenitor. Starting from a 
host halo at 2 = 20, we iterate this procedure back in time, 
thus obtaining a complete merger tree down to the mass 
resolution fixed at 10 particles per halo. 

Fig. 1 shows a schematic representation of a merger tree 
to help us define the terminology used throughout the paper. 
The main branch of the merger tree is defined as the branch 
tracing the main progenitor of the main progenitor of the 
main progenitor etc... (i.e. the branch connecting the dark 
grey haloes in Fig. 1). We will use the term satellites to refer 
to all progenitor haloes accreted by the main progenitor, 
and donating at least 50% of their mass to the host halo 
at 2 = zq. In Fig. 1 these are the light grey haloes A, B, 
C and D (as defined before they were accreted onto the 
main branch) . Extending the same definition to each of these 
satellites haloes, we can go back in time and trace for each 
satellite its own main branch and satellites. For example, 
satellite D has accreted its own satellites a, b and c at some 
redshift 23, Z2 and zi, respectively. 

We stress that not all satellites will survive to 2 = zq: 
many of them will be destroyed by tidal effects, encounters 
and other dynamical processes. As an example, satellite a 
in Fig. 1 only survives until 2 = 22. For this reason, in order 
to identify substructures at zq, at any level of the hierarchy, 
we need to climb the merger trees of all satellites (at any 
level of the hierarchy) only up to satellites that will make 
it to 2 = 20, but not further. In the example of Fig. 1, 
we will stop climbing the merger tree of satellite D at 23, 
which is the lowest redshift when any possible substructure 
present inside D at 2 = 23 has no self-bound counterpart 
at 20 . We will call D(zz) a heirless satellites. In the same 
spirit, the other heirless satellites in Fig. 1 are A(z4), B(z^), 
C(z2), b(zz) and c(z2). Their self-bound counterparts at 20 
constitute the substructure population of the host halo, at 
all levels of the hierarchy. In what follows we use subhaloes'^ 
to refer to the i th level of this substructure hierarchy, with 
i — 1 corresponding to the first order of subhaloes, i.e. the 
ones accreted directly by the main progenitor. Hence, in 
Fig. haloes A, B and C are subhaloes^ 1 ', while b and c 
make up the population of subhaloes^ 2 ' of the host halo at 
z — zq. We will use the term 'substructures' to refer to 
subhaloes of the entire hierarchy (i = 1,2,3, ...). The mass 
of a substructure at 20 is defined by the clump of particles, 
originally belonging to a heirless satellite, that at 2 = 20 
are still self-bound within the tidal radius of that clump. 



This procedure is iterative, because the tidal radius itself is 
defined using the self-bound mass of the clump. 

To implement this algorithm in our code, assume that 
we have a simulation with N snapshots, and let Zi be the 
redshift of snapshot i. Let 2jv and 20 be the redshifts of 
the earliest and final snapshot, respectively. Starting from 
the population of subhaloes identified by SURV at 20, we 
proceed along the following steps: 

• given a subhalo s(zq), accreted as satellite at Zi > 20, 
we climb one step along its own merging history tree and 
read the information about its progenitors at Zt+i; 

• for each of these progenitors, we trace the position of its 
particles to 20, and see if there is a self-bound clump within 
its own tidal radius (as defined in Tormen et al. 1998); 

• if at least two progenitors identified at m+i (the main 
progenitor and another satellite) have a counterpart clump 
s'(zo) made of at least 10 self-bound particles, this means 
that 5(20) has substructures inside itself. We therefore up- 
date the previous subhalo catalogue by replacing s(zo) with 
the counterparts s'i(zq), s' 2 (zq), of the survived progeni- 
tors. Notice that one of these progenitors will always be the 
main progenitor of s(zo), so this procedure will certainly re- 
identify the core of s(zo). Any other counterpart s'i(zo) will 
identify a new subhalo at a further level of the substructure 
hierarchy. 

• if, of all progenitors identified at z»+i, only the main 
progenitor has a self-bound counterpart at zq, this means 
that s(zo) has no further levels of substructure, and the 
(satellite) main progenitor of s(zq) at Zi is a heirless satellite. 
In this case we leave the subhalo s(zq) as is, and proceed to 
the next one in the catalogue. 

After we have scrutinized all subhaloes once, we start 
again and consider the updated subhalo catalogue (which 
now contains new subhaloes at one further level of hierar- 
chy). We repeat the steps described above, skipping however 
the subhaloes coming from heirless satellites, and split any 
new subhalo as described. We iterate the whole procedure: 
each time we replace the subhaloes - which at zq may be 
split in two or more self-bound sub-units - with the sub- 
units themselves. At some point, all subhaloes at zq will be 
the counterpart of some heirless satellite. We call the result- 
ing population a catalogue of substructures. This catalogue 
is the one we use for our investigation. 
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Figure 2. Dark matter particle distributions for the most massive halo (AT = 1.85 x lO^fe^M©) in the GIF2 run at z = 0, the radius 
of the sphere is 2.54/i _1 Mpc. Top left: all dark matter particles within the virial radius. Top right: diffuse dark matter not associated 
with any substructure. Bottom left: substructures, i.e. self-bound particles belonging to heirless satellites haloes (see the main text for 
more details). The particle distribution in the top right and bottom left sum to give that in the top left. Bottom right: subhaloes' 1 ) , i.e. 
self-bound particles from satellite haloes accreted directly by the main progenitor (Giocoli ct al. 2008a). The two bottom panels only 
show substructures and subhaloes' 1 ) with at least 10 particles. Clumps at radii r < 0.05iJ v i r are not well denned, and were excluded 
from the study by Giocoli et al. (2008a) and hereafter. 



4 RESULTS 

Figure shows the dark matter distribution for the most 
massive halo found in the GIF2 run at z = 0. The top 
left panel shows all particles inside the virial radius. The 
top right panel shows halo particles not bound to substruc- 
tures (dust particles). The bottom left panel shows particles 
bound to substructures ( i.e. the final result of our proce- 
dure). The particle distribution shown in the top left panel 
is the sum of that in this panel and that shown in the top 
right. The bottom right panel shows the particles bound to 
subhaloes' 1 ' (i.e. the starting point of our procedure); these 



are the objects identified by Giocoli et al. (2008a). A com- 
parison between the two lower panels shows that the largest 
subhaloes' 1 ' are split into smaller and smaller substructures 
as we proceed further along the merger tree, until we reach 
the heirless satellites. Note that the total mass in the two 
bottom panels is different for two reasons. First, some of 
the particles in the Giocoli et al. (2008a) subhaloes are no 
longer bound to any of the substructures. Second, some sub- 
structures found by our new algorithm (bottom left) were 
not identified as subhaloes' 1 ' (bottom right). While the old 
algorithm SURV only follows the most massive piece, our 
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new algorithm identifies all these pieces as satellites of the 
original system. Figure Al shows the net change in the mass 
fraction associated with substructures. 

In what follows we discuss the substructure of the host 
haloes as function of mass and redshift. We consider five 
values of zq: 0, 0.5, 1, 2 and 4, and for each, we consider 
seven mass bins. Throughout we only use host haloes whose 
main progenitor never exceeds the final mass M ZQ by more 
than 10%. Table lists the number of haloes in each of 
these bins. Note that substructures in the innermost region 
(within 5% of the virial radius) are excluded from further 
study, as their identification is uncertain due to the extreme 
density of the host. 

The resulting substructure mass functions are shown 
in Figure 3; masses are in units of the virial mass of the 
host halo. The various symbols and line types in each panel 
show results for different host masses. The solid curve (same 
in each panel) shows the unevolved subhalo mass function 
(equation 2 of Giocoli et al. 2008a), i.e., the mass function 
of progenitor haloes accreted onto the main branch of the 
merger tree: 



-2 



dN(m sh \M) 
d In m s b 



No(M zo )C exp(-/3£ 3 ) =£/(£): 
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with ^ = m sb /M Z0 , a = -0.8, f) = 12.2715 » 0.43" 3 and 
normalization No(M Zo ) = 0.18. Note that the normalization 
is independent of both the mass of the host halo and the 
redshift at which it was identified even though our notation 
suggests otherwise. The reason for our notation will become 
clear shortly. For comparison, the dashed line shows a power 
law with slope a = —0.9 and normalization No = 0.03: it is 
steeper than the solid line. Performing least squares interpo- 
lation on the linear trend of the substructure mass functions 
we notice that this value of the slope is in agreement with 
the value estimated in the more massive bins, while it tends 
to overestimate a little the slope in the smaller ones. This 
because subhaloes in the latter are mainly preserved and not 
split in sub-subhaloes. 

Giocoli et al. (2008a) have shown that, even after tidal 
stripping, the evolved mass function of subhaloes' 1 ' (i.e. for 
the objects in the bottom right panel of Figure 2) has the 
same shape as the unevolved one (e.g., power law slope 
a = —0.8 at small m s b/M Z0 ). I.e., the mass function for 
objects like those in the bottom right panel of Figure 2 is 
given by sliding the solid curve downwards by an amount 
that depends on M Z() (and zo). The reason for this is not 
completely straightforward, since mass-loss alone would be 
associated with sliding to the left, not down. As discussed 
in van den Bosch et al. (2005), this happens because the 
mass loss depends on the time spent by a satellite within the 
potential well of its host, and Giocoli et al. (2008a) showed 
that the mass function of objects which accreted before and 
after the main progenitor had acquired half the final mass of 
the halo are both described by equation (2), but with differ- 
ent normalization factors (lower by factors of 0.57 and 0.43 
respectively). Crudely speaking, all the subclumps accreted 
prior to the half-mass time have been erased completely, and 
all those accreted later are still present, but they have been 
stripped. So, crudely speaking the evolved mass function 
should be given by taking equation ( ), but with normal- 
ization lower by a factor of 0.43, and then shirting it to the 
left. To match amplitudes, the typical shift is about a factor 



O 
X 



o 



-3 - 



-4 



-5 - 



j — i — i — i — r - 

fit a= -0. 

mean N. 




11 12 13 

!og(nn sb ) 



Figure 4. Substructure mass function per unit host halo mass 
(in units of 10 10 h~ 1 Mq). Different data points refer to the same 
mass bins as in Figure 3. The solid line shows a power law distri- 
bution with slope a = —0.9 and a mean normalization computed 
over the different mass bins (see the text for more details) . The 
dashed lines show the fit to the data for each mass bin plus the 
high mass exponential cut-off. 



of 3. However, because low mass objects typically formed 
at higher redshifts, their subhaloes have suffered more mass 
loss, so the shift is larger for smaller M ZQ . Over the mass 
range where the mass function resembles a power law, the 
leftwards shift resembles a downwards shift in the normaliza- 
tion of equation (2), so the number of subhaloes of a given 
m B b/M Z0 is larger for larger M ZQ . In the next section, we 
show that the scatter around this relation is also directly 
related to the halo assembly process. 

The symbols in Figure 3 do not show this evolved sub- 
halo mass function. Rather, they show the substructure mass 
function associated with our new algorithm (i.e., for the ob- 
jects shown in the bottom left panel of Figure ). Notice that 
this function has a steeper slope (a = —0.9) than that for 
the subhaloes' 1 ' (at least in the most massive host haloes). 
However, equation (2) still provides a good description of 
the full shape, if we simply set a — —0.9, keep j3 — 12.2715, 
but now allow the normalization No(M ZQ ) to depend on the 
mass (and, we show later, on zq) of the host halo. 

Figure I shows the substructure mass function if one 
does not normalize all m s b by M ZQ . Namely, it shows 



AN M 



fn.sb 



AN 



dlnm sb M dm sh 
where the normalization factor is 



= Nm m"h exp (-Pf) , (3) 



N 



M 



No(M zo ) 
M 1+a 



(4) 



The dashed curves show this expression with a — —0.9 and 
normalization determined by chi-square minimization to the 
data for each bin in Mo- Figure 5 shows that, with the ex- 
ception of the lowest mass bin, log Nm — —3.03 almost 
independent of mass. We believe the lowest mass bin, for 



The Substructure Hierarchy in Dark Matter Haloes 7 




log(m s5 /M Zo ) 

Figure 3. Substructure mass functions at five different redshifts zq as indicated. We divide the host haloes in seven different mass bins 
and plot the mass function (in units of the host halo mass M ZQ ) at that redshift. In each panel the dashed line shows a power law with 
slope a = —0.9 and normalization No = 0.03 that fits well the data points at zg = 0. For comparison, the solid curve represents the 
unevolved subhalo mass function (equation 2). 



which the haloes have the fewest particles, is compromised 
by discreteness effects: the substructure mass loss rates of 
the objects in this bin are rather discontinuous. 



The substructure population in this GIF2 simulation 
has also been studied by Gao et al. (2004). They used 
a different post-processing pipeline and their algorithm to 
identify substructures, SUBFIND, is different from ours. 
So it is remarkable that our mass functions agree rather 
well: we both find a = —0.9. However, our normalizations 
are different: they find —3.2 whereas we find —3.03. This 
can be traced to the fact that we define haloes as being 
A v i r (z)p cr j t (z), whereas they fix A v i r = 200. In effect, their 
haloes are smaller and less massive than ours (Appendix 
shows that the volume V200 ^ 0.4 Kdr at z = 0). Hence, at 
the same numerical value for the halo mass, our haloes are 
'effectively' more massive than theirs, which means smaller 
formation redshifts, and hence higher mass fractions of sur- 
viving substructure. We believe this to be responsible for 
the difference in normalization. Appendix A gives a more 
detailed discussion of the differences between our new algo- 
rithm, that of Gao et al. (2004), and that of Giocoli et al. 
(2008a). 



CP 
O 




log(M ) 



Figure 5. Host halo mass dependence of the normalization factor 
iVjvf in eq. (3) at zg = 0. The error bars were obtained by boot- 
strapping the host halo sample in each mass bin and computing 
the rms of the normalization factors obtained for the bootstraps. 
The dashed line in the top panel shows the mean weighted by the 
errors of the data points excluding the smallest mass bin (see the 
main text for more details). 
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5 MASS FRACTION AND SUBSTRUCTURES 

In this section, we show how the substructure mass fraction 
correlates with other properties of the halo (mass, concen- 
tration, formation time). Most of these are expected - our 
results quantify the trends. We then study how the number 
of substructures depends on the redshift at which the host 
halo is identified. 



5.1 Correlations with mass, formation time, and 
concentration 

The structure of a dark matter halo is related to its accre- 
tion history (Zhao et al. 2009, and references therein). To 
quantify the different halo assembly histories we will use the 
formation time, defined as the earliest time when the main 
progenitor reaches a mass that is half of the final mass of the 
halo (e.g., Lacey & Cole 1993). Previous work has shown 
that massive objects formed more recently (e.g., Lacey & 
Cole 1993; van den Bosch 2002; Sheth & Tormen 2004; 
Giocoli et al. 2007). Systems that form at high redshift, 
when the universe was denser, are denser. This, with the 
fact that massive objects formed more recently gives rise to 
the mass concentration relation (Navarro et al. 1996; Bul- 
lock et al. 2001; Maccio et al. 2007; Neto et al. 2007; 
Zhao et al. 2009) and its scatter. In addition, using ex- 
tended Press-Schechter theory, van den Bosch et al. (2005) 
showed that the unevolved subhalo mass function is uni- 
versal, which was confirmed with numerical simulations by 
Giocoli et al. (2008a). Hence, they predicted that less mas- 
sive haloes, which on average form earlier (i.e., accrete their 
subhaloes earlier), would have less surviving substructure at 
the present, simply because their substructure would have 
been exposed to tidal stripping for a longer duration, van 
den Bosch et al. (2005) thus predicted a correlation be- 
tween f a and mass, but also that the scatter in /„ at fixed 
mass be correlated with the formation time of the halo. 
To check this, for each halo we computed 



and, for each mass bin, 



Mo 



(fs) 



Nhaloe 



When computing these sums, we considered all substruc- 
tures with at least 10 self-bound dark matter particles that 
were more than 0.05i? v ir from the center of the host halo. 

Figure 6 shows this mean substructure mass fraction at 
z — for a range of host halo masses. The open circles show 
results for the full set of haloes in each mass bin, whereas 
triangles and diamonds show the subset which formed below 
and above the mean formation redshift (for that mass). At 
z = 0, we have 13467 (7158), 4407 (3898), 1834 (1515), 610 
(576), 215 (246), 72 (55), 18 (17) and 2 (2) haloes with a 
formation redshift higher (lower) than the mean, for our 
seven mass bins. This shows that the normalization factor 
No of the substructure mass function depends on Mo and on 
the formation redshift zj. For a given mass, lower formation 
redshift translates into a higher substructure mass fraction 
and smaller smooth dark matter component Mo(l — / s ). The 
three dashed lines show least squares fits to the three sets 
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Figure 6. Fraction of present-day host halo mass in substruc- 
tures. Open circles show the substructure mass fraction consid- 
ering all haloes in each of the seven mass bins. Filled triangles 
and open diamonds represent the mean substructure mass frac- 
tion considering only haloes with a formation redshift lower and 
higher than the mean in the corresponding mass bin, respectively. 
The error bar shows the rms. The dashed lines show least squares 
fits to the three sets of data points. 



of data points: 



(f s )=alog(M ) + b, 



(5) 



where a and b for the lines going from the top to bottom in 
Figure 6 we have: 

f (0.020 ± 0.004, -0.16 ± 0.06) if z f < z f , 
(a, b) = < (0.024 ± 0.001, -0.24 ± 0.02) for all haloes, 
[(0.027 ± 0.002, -0.32 ±0.03) if z f > z f . 

Figure shows another way of presenting the joint 
distribution of halo mass, formation time and f 3 . In this 
case, we study the correlation between f s and Zf for differ- 
ent masses. The black diamonds show this correlation for 
all haloes more massive than 10 Mg/A, the dashed line 
shows a least squares fit: 

log(/ s ) = (-3.133 ± 0.13) log(l + z f ) - (0.48 ± 0.05), (6) 

and the other symbols show the correlation for narrow mass 
bins. Notice that the substructure mass fraction is tightly 
correlated with the halo assembly redshift for haloes of all 
masses, indicating that, to good approximation, we can write 
that f s — f 3 (zf). The fact that this f„ — Zf correlation de- 
pends little on mass Mo, whereas the f s — Mo correlation 
depends strongly on Zf indicates that the f s — Mo correla- 
tion is almost entirely determined by the fs — Zf and zj — Mo 
correlations. From the figure we notice also that for haloes 
formed quite recently, a high value of the substructure mass 
fraction indicate that they have undergone a recent almost 
1:1 major merging event. This can also been understood if 
we consider the conditional distribution of masses at forma- 
tion studied by Sheth & Tormen (2004). In this work they 
show that the distribution is tightly peaked around 1/2 for 
large formation redshift (since there were few progenitors of 
large mass at high-z), whereas it is much broader for lower 
zf. so low Zf means more recent 1:1 mergers. 
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Figure 7. Correlation between the substructure mass fraction 
and host halo formation redshift. Different symbols show results 
for the mass bins as in Figure . Filled diamonds show the relation 
averaged over all haloes more massive than 10 n - 5 M Q /h at z = 0; 
dashed line shows a least squares fit. 



Figure 9. Concentration— formation redshift correlation. Differ- 
ent symbols and line styles show results for different mass bins (as 
in Figure 7). Dashed line shows a least squares fit to the relation 
traced by the filled diamonds, which show the average over all 
haloes at z = 0. 
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Figure 8. Mass-concentration relation at z = Q. Open circles 
show the mean for the full sample at each mass; filled triangles 
and open diamonds show this relation for haloes with formation 
redshifts that are lower and higher than the mean (similar to 
Figure ). Dashed lines show least squares fit to these relations. 
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Figure 10. Dependence of mass fraction in substructure on host 
halo mass and concentration. Open circles show the full sample, 
triangles and stars show the average over systems with above and 
below average concentrations for the mass bin. Dashed lines show 
least squares fits to these trends. 



Halo formation times are not observable, so it is in- 
teresting to see if halo substructure correlates with other 
structural properties. Halo concentrations are potentially 
observable, and are known to correlate with mass and for- 
mation time (e.g. Navarro et al. 1997; Wechsler et al. 2002; 
Zhao et al. 2009). Figures 8 and 9 show these relations in 
our simulation. More massive haloes are less concentrated; 
at fixed mass, the objects which formed at higher redshift 



are more concentrated. Different data points refer to various 
mass bins, as in previous figures, and error bars represent the 
rms around the mean value. We estimate the halo concentra- 
tion Cvir = R v ir/r s fitting the spherical density distribution 
with an NFW profile. In Figure 8 the three lines show 

log( log(Mo) + b (7) 
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Figure 11. Dependence of mass fraction in substructure on host 
halo concentration and mass. Symbols and line styles are for dif- 
ferent host halo masses (same as Figure 7). 

with 

f (-0.102 ± 0.008, 2.26 ± 0.10) if z f < %, 

(a, 6) = I (-0.114 ± 0.006, 2.49 ± 0.08) for all haloes, 

( (-0.122 ± 0.005, 2.65 ± 0.07) if zj > z f . 

In Figure 9 the dashed line shows the least squares fit to 
all masses. However because the normalization of this cor- 
relation depends only weakly on host halo mass, it fits the 
smallest mass bin, which dominates the numbers. 

Figure shows the result of remaking Figure 6, but 
now replacing host halo formation time with concentration. 
This shows that, at fixed halo mass, more concentrated 
haloes have smaller substructure mass fractions. The three 
dashed lines show the least squares fits to the data points 
where 

{(0.251 ± 0.005, -0.24 ± 0.07) if c vir < c vir , 
(0.024 ± 0.001, -0.24 ± 0.02) for all haloes, 
(0.026 ± 0.003, -0.29 ± 0.03) if c vir > c vir . 

Figure shows a similar remake of Figure 7; f s de- 
creases as concentration increases, approximately indepen- 
dent of halo mass. If we ignore the smallest mass bin, which 
is most affected by the mass resolution of the simulation, 
then data points are well-fitted by an inverse correlation be- 
tween the mass fraction and the concentration of the host 
halo. 

5.2 Redshift dependence 

At fixed host halo mass the number of substructures depends 
on the redshift at which the halo was identified: it is larger 
at higher redshift. Although they did not emphasize this 
aspect, Angulo et al. (2008) showed that the subhaloes of 
hosts identified as being 200 times the critical density in 
the Millennium simulation using SUBFIND show a similar 
trend. 

Substructure abundance depends on how much strip- 
ping has occurred, since this is what changes the univer- 
sality of the unevolved subhalo population. This stripping 



depends on how long a satellite has spent inside its host, 
and on the dynamical timescale within the host (van den 
Bosch et al. 2005; Giocoli et al. 2008a). These timescales, 
that represent the time over which changes in one part of a 
body can be communicated to the rest of that body, depend 
on when the host halo is identified (rd yn oc p v lJ 2 ), so this 
may be why the higher redshift hosts have more substruc- 
tures. If we use T to denote the ratio of these two timescales 
(T = (t(z) — £ a cc) /idyn(z)), then large values of T should be 
associated with more stripping. The top panels of Figure 12 
show the distribution of T for substructures in different host 
halo masses (right to left) identified at two redshifts, zo — 
and 2 (solid and dashed histograms). This shows that, for 
host haloes observed at higher redshift the distribution of T 
peaks at lower values. Substructures in these systems that 
have spent less time in the potential well of their host suf- 
fer less tidal stripping than those in present-day hosts (of 
the same mass). The bottom panels show the correlation 
between the survived mass fraction m B \ D (zo) /m v i T (z sxcc ) and 
T. As discussed in van den Bosch et al. (2005) and Giocoli 
et al. (2008a), when the instantaneous mass loss rate does 
not depend on the host halo mass, then the retained mass 
fraction decreases as exp (— T). 

5.3 Poisson model for substructure abundance 

Figure shows how the mean number of substructures 
scales with host halo mass, again only counting substruc- 
ture with ten particles or more (i.e., m s b. m i n = 1-73 x 
10 10 7i _1 Mq). The dashed line shows the result of fitting 

(N a ) = A)(m s b, mm) Mq , (8) 

to the measurements. A least squares fit returns /3 = 0.97 ± 
0.03 and log(Aj) = -11.79 ± 0.34. Of course, this normal- 
ization factor Aq depends on the smallest substructure mass 
considered, but /3 « 1 is consistent with the fact that the 
normalization constant Nm m equation (3) is independent 
of Mo. The other panels of Figure L3 show the second and 
third factorial moments of this distribution: {N S (N S — 1)) 
and {N S (N S — 1)(N S — 2)). If the scatter around the mean 
were Poisson, these would equal {N s } 2 and (N s ) 3 respec- 
tively. Evidently, our substructure counts deviate from the 
Poisson model only slightly (by a factor of ~ 5% at large 
masses, and slightly increasing toward small ones). Inter- 
estingly, a similar trend was found by van den Bosch et al. 
(2005) using merger trees constructed using extended Press- 
Schechter theory. 

Figure shows how the mean number of satellites 
(more massive than 1.73 x 10 10 /i _1 Mg) depends on the red- 
shift at which the host halo was identified. Different symbols 
show counts in haloes at zo = 0, 0.5, 1, 2, 4. Notice that, for 
the reasons discussed earlier, at a given host mass, the mean 
counts are higher for the hosts identified at higher redshift, 
but that the curves are parallel to one another. The higher 
number of substructures at high redshift can also be related 
to an increase in the galaxy merger probability (Wetzel et al. 
2009b; Hester & Tasitsiomi 2009) and have implications for 
galaxy formations and mergers. We showed previously that 
the substructure mass fraction f s is anti-correlated with halo 
concentration c, approximately independent of halo mass. 
For a given host halo mass, c is smaller at high redshift (e.g. 
Zhao et al. 2009). Hence, if the f s — c anti-correlation is the 
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Figure 12. Top panels: Distribution of the time spent by substructures in their host haloes (i.e. since they were accreted), in units of 
the dynamical time scale at zq = (solid) and zq = 2 (dashed). The three panels show results for three bins in host halo mass. Bottom 
panels: Retained mass fraction for substructures which have survived within their host haloes as a function of (scaled) time, for the same 
two choices of zq an d for the same three bins in host mass. 
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Figure 13. First three factorial moments of the substructure counts as a function of host halo mass. Only substructures with more than 
10 self-bound particles, and further than 0.05i? v ir of the host halo center were counted. A Poisson distribution would have value unity 
in the middle and right-most panels. 
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Figure 14. Dependence of the mean number of satellites more 
massive than 1.73 X 10 10 /i _1 Mq on halo mass, and on the red- 
shift at which the halo was identified. At fixed mass, host haloes 
identified at higher redshift have more substructures. 




1 +z 

Figure 15. Redshift evolution of the concentration of an M*- 
halo (here shown in units of the concentration of M, haloes at 
2 = 0) as predicted by the model by Zhao ct al. (2009). 

same whatever the redshift at which the host halo is iden- 
tified, then one expects f s to be larger at high redshifts, in 
qualitative agreement with the trend shown in Figure 14. 

The hypothesis that the f s — c anti-correlation does not 
evolve allows us to be more quantitative, provided that we 
have a model for how halo concentrations depend on the 
redshift at which they were identified. Zhao et al. (2009) 
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Figure 16. Same as previous figure, but now rescaled by a factor 
which reflects the anti-correlation between substructure fraction 
and halo concentration, and the redshift dependence of the halo 
concentration. 

describe an accurate model for this, in which a halo's con- 
centration is related to the time at which its main progenitor 
first assembled 4% of its final mass. We used their model to 
make Figure 15, which shows the concentration of an M 4 
halo as a function of redshift z. (I.e., M, represents the typ- 
ical halo mass at redshift z, defined by S{M*) = 5 c {z) 2 , 
where S c (z) is the critical overdensity required for spherical 
collapse at z and S(M) is the variance in the linear fluctua- 
tion field when smoothed with a top-hat filter which contains 
mass M. Thus, M* is smaller at higher redshift.) This shows 
that, between z = and 4, c(M« (z))/c(M, (0)) oc (l+z)~ 1/2 : 
although M,(z) is smaller at high redshift, so too is its con- 
centration. 

Since f s oc c -1 at z — (Figure 11), the hypothe- 
sis that the f a — c anti-correlation does not evolve suggests 
that the typical value of f a at some zq will be proportional 
to c(M t (zo))~ 1 ■ Therefore, subhalo counts at Zq should lie 
above those at z — by a factor of [c(Mq)/c(M* )] ~ 
(1 + zq) 1 / 2 . Figure 16 shows the result of dividing the zo 
counts shown in Figure 14 by (1 + zq) 1 ^ 2 : this does indeed 
bring the counts at different zq into good agreement with 
one another. To highlight how well this works, Figure 17 
shows the analogue of Figure 5 but for these rescaled counts 
(symbols and line types as in Figure 14): the rescaled counts 
are much less redshift-dependent, suggesting that our crude 
accounting for the effects of the f s — c correlation on subhalo 
abundances is reasonably accurate. 



6 GLOBAL MASS FUNCTION OF 
SUBSTRUCTURES 

The mass function of substructures, integrated over all host 
haloes, plays an important role in models which relate sub- 



The Substructure Hierarchy in Dark Matter Haloes 13 



-2.9 - 



<? -3.1 



-3.2 




1 2 



13 

log(Mj 



1 4 



15 



Figure 17. Rescaled normalization factor for the number of sub- 
structures per host halo mass, shown in the same format as Fig- 
ure 5. Different symbols are for host haloes identified at different 
redshifts zq (as in Figure 14). 
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structures to galaxies (Cooray & Sheth 2002) and to the ex- 
pected gravitational lensing signal on small scales (Sheth & 
Jain 2003; Maccio & Miranda 2006; Natarajan et al. 2007; 
Meneghetti et al. 2007). The global substructure mass func- 
tion can also be useful for modeling the 7-ray background 
due to dark matter particles annihilation (Fornengo et al. 
2004; Giocoli et al. 2008b, 2009). We now have the necessary 
ingredients to estimate this quantity: 



dn(m sb ) 
d In m sb 



dN(m sb \M) An(M) 
d In m sb AM 

= N Mo (l + z) 1/2 (m sb ) a 

Ja , dn(M) , 
AM — K - ' exp (- 
d In M v 



AM 



(9) 



where An /AM denotes the comoving number density of dark 
matter haloes. The histograms in Figure 18 show the re- 
sult of performing this integral by using equation (2) with 
a — —0.9, /3 = 12.2715, the rescaled normalization from Fig- 
ure 17 for AN(m\M)/Am, and the Sheth & Tormen (1999) 
expression for the mass function of the hosts, An/AM, at 
z — 0, 0.5, 1, 2, 4 (top to bottom). The short-dashed curves 
show An/AM at these same redshifts. 

Because of the exponential cutoff in the subhalo 
mass function, equation (9) cannot be integrated an- 
alytically. However, we have found that j/(m sb ) = 
m sb dn(m s b)/dlnm s b is well fit by 



y(m s h) = A m^ b exp 



m sb 
m 



(10) 



A similar functional form has been used to fit galaxy velocity 
dispersions (Sheth et al. 2003), so our parametrization is 
intended to help clarify the connection between galaxies and 
halo substructure. To determine the parameters Ao, n, 7 and 
mo, we perform a least-squares fit to the log of the counts 
at the small mass end: this determines the slope r\ and zero- 
point ao of the power law at small masses. We set Ao = 10 a ° , 
and we estimate the other two parameters, 7 and mo, by 



Figure 18. Number density per M 2 /p of substructures (his- 
tograms) at redshifts 2 = 0, 0.5, 1, 2, 4 (top to bottom at massive 
end); dashed curves show the result of fitting to equation (10). For 
comparison, short-dashed curves show the host halo mass func- 
tion (from Sheth & Tormen 1999) at the same redshifts. The 
solid line show the halo+substructure mass function and the dot- 
dashed the analytical solution for the global substructure mass 
function when /5 fs (eq. 12). 



minimizing 



log(j/(m s b)) - log m sb 



dn(m sb ) 
' d In m sb 



(11) 



where the sum is over the various values of In m sb at which 
we computed the integral (the centers of the histograms 
shown in Figure 18). Table 2 lists the resulting best fit pa- 
rameters. As can be seen from the last column of the Table, 
the shape is fit to an accuracy of about 2%. Solid curves in 
the figure show the halo+substructures mass function. 
0, equation (9) becomes 

A 



When j3 : 
dn(m sb ) 



d In m sb 



N Mo (l + z) 1/2 (m sh ) a p 



x [0.812r(0.2,0.353^) + T(0.5, 0.353^)] (12) 



where A w 0.322 is the normalization of the halo mass 
function and v = 5 c (z) 2 / S(m). The dot-dashed curves in 
Figure show this approximation. Clearly, neglecting the 
high-mass exponential cut off is not a good approximation. 



7 SUMMARY AND CONCLUSION 

We estimated the mass function of substructures by follow- 
ing all the branches (rather than just the main branch) of the 
merger history tree of a halo. This distribution, diV/dm sb , 
depends on the mass of the host halo, when the halo was 
identified, and its formation time. However when the num- 
ber of substructures per host halo mass is considered, the 
distribution turns to be "universal" at a given redshift, but 
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Table 2. Best fit parameter to the global substructure mass function at the five considered rcdshifts. The last column shows the rms 
scatter between the distribution and the best fit. 



simulations with higher mass resolution are required to de- 
termine if this is true even when the host halo mass is smaller 
than 10 11 ft -1 Mq. This universal shape is characterized by a 
power law of slope —0.9 at low masses, times an exponential 
cut-off when the substructure mass approaches that of the 
host halo. This distribution can be integrated to obtain the 
substructure virial mass fraction and the mean number of 
clumps at a given host halo mass. 

At fixed host halo mass, the haloes which formed at 
higher redshift tend to have less substructure. In addition, 
more massive host haloes have more substructures than less 
massive ones. The first trend arises because satellites which 
spend more time in the potential well of their host lose a 
larger fraction of their initial mass (or can be totally dis- 
rupted). Massive haloes formed more recently, so this same 
mechanism explains the second trend. Similar trends are 
seen if we replace halo formation time (which is not observ- 
able) with the concentration parameter of the halo density 
profile (which is). 

We also showed that the Poisson model, with a mean 
that depends on host halo mass, is a reasonable descrip- 
tion of the substructure counts within a host halo. However, 
the second and third factorial moments of the substructure 
counts are ~ 5 percent larger than for a Poisson model with 
the same mean counts. As observations improve, this may 
become important in Halo Model parametrization of large 
scale structure which relate halo substructure to galaxies. 

At fixed mass, the mean number of substructures above 
some minimum mass is larger at high redshift. Although this 
is not unexpected - the higher redshift haloes have had less 
time to destroy their substructures - we argued that the red- 
shift dependence could be well-approximated by assuming 
that the anti-correlation between substructure abundance 
and host halo concentration we measure at z = does not 
evolve. At fixed mass, higher redshift haloes are less concen- 
trated - in effect, our assumption is that the same physics 
which affects the concentrations also affects the substruc- 
ture. This allows us to provide a simple prescription for the 
redshift evolution of substructures, which we use to provide 
analytic approximations for the result of integrating the sub- 
structure mass function over all host halo masses. These re- 
sults should find use in halo model interpretations of the 
(evolution of the) galaxy clustering signal in large scale sky 
surveys (e.g., van den Bosch et al. 2007; Wake et al. 2008). 
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APPENDIX A: COMPARISON WITH 
PREVIOUS WORK 

Our new substructure algorithm is a modification of that in 
Giocoli et al. (2008a) and is rather different from that of 
Gaoetal. (2004). 

Our new algorithm returns smaller substructure mass 
fractions, especially at large host halo masses, than did our 
old one (i.e., that of Giocoli et al. 2008a)); it is more conser- 
vative about identifying self-bound structures. This is shown 
in Figure Al, where the change at the high mass end is more 
than 50%. It is smaller at the low mass end, where the host 
halo masses are approaching that of the mass resolution of 
the subhaloes. This shows explicitly that if one wants to 
think of our new algorithm as taking the subhaloes of Gio- 
coli et al. (2008a) and partitioning them up into smaller 
substructures, then some of the subhalo mass from the old 
algorithm is now assigned to the host halo rather than to 
substructures. The new value of the substructure mass frac- 
tion for the most massive haloes agrees pretty well with 
other algorithms, that give values around 10% of the host 
halo mass (Springel et al. 2001; Weller et al. 2005). 

The changes in mass fraction between the two algo- 
rithms shown in Figure Al can be understood also consid- 
ering the fact that when we split a subhalo in pieces (sub- 
subhaloes), we must also recompute the centre-of- velocity 
(from which we compute the kinetic energy) and, more im- 
portantly, we calculate the new potential energy only using 
the new subsets of particles. This forces the new potential 
energy to be less negative than before (i.e. the same particle 
is assigned a less negative potential energy than before, be- 
cause the potential wells are shallower), and so many parti- 
cles which reside in-between the density peaks of the various 
sub-subhaloes, even if they were bound to the subhalo as a 
whole, happen to be unbound to any of the sub-subhaloes, 
and often are unbound to the main subhalo itself. This is be- 
cause the splitting removed all the secondary density peaks, 
so the global potential well becomes shallower decreasing 
the self-bound mass of the clumps. 

Figure A2 compares the substructure mass functions 
found by our algorithm, that of Gao et al. (2004), and the 
subhalo mass function of Giocoli et al. (2008a). To facilitate 
comparison with Gao et al. (2004), we have chosen to use 
the same mass bins and host halo definition as Gao et al. 

1 Comparing our new algorithm with the result of SUBFIND for 
the most massive halo at z = we found respectively 6.5% and 
6.2% of A^oo.c to be in substructures. 



0.15 
v 0.1 



0.05 







— A subhalos 

H O substructures 


I 


1 1 






i ~ 


















f 


r" 




ri 




rill 


. 1 






i i 


12 


13 


14 






15 



log(M ) 

Figure Al. Mass fraction in the subhaloes of Giocoli et al. 
(2008a) and in substructures found by our new algorithm (filled 
and open symbols, respectively). 
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Figure A2. Comparison of substructure mass functions, ex- 
pressed as a function of rn s ^/M200 found by our algorithm, that 
of Gao et al. (2004) (SUBFIND), and the subhaloes' 1 ) of Giocoli 
et al. (2008a). 

(2004); i.e., host haloes are 200 times the critical density. We 
will use M200 do distinguish this definition from Af v i r , the 
quantity used in the main text. Note that M200 < M V - 1T , with 
the difference increasing at late times (lower redshifts). Our 
new algorithm (open symbols) takes the subhaloes returned 
by our old algorithm (triangles) and partitions them up into 
smaller pieces; it removes objects from large m s b ~ M200 
and adds objects to m s b <C M200. The resulting distribution 
is indeed very similar to that found by Gao et al. (2004) 
using SUBFIND (crosses). The dashed and the dot-dashed 
lines show the fit-by-eye of the cluster size halo by De Lucia 
et al. (2004) and Gao et al. (2004). 
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Figure A3. Substructure mass functions found by our algorithm 
when host haloes arc defined to be 200 times the critical density. 
This is the analog of Figure 1 in the main text. Solid line show 
equation (1) by Gao et al. (2004). 



Figure , which shows diV/dln(m s b /M200), is the ana- 
log of Figure '■ in the main text. However, in our discussion 
of (diV/dmm s b)/M v j r in the main text, we found that the 
normalization of our substructure mass function appeared 
to be larger that that reported by Gao et al. (2004), and 
we argued that this could be attributed to the difference in 
how we define host halo masses. Since the definitions are the 
same here, Figure A3 shows (diV/dlnm s b)/M2oo, which is 
the analog of Figure I . Note that the normalization is indeed 
lower for M200 than it is for M v i r . 

In addition, we argued in the main text that because 
M200 < M v i r , the host halo formation times for haloes of a 
given M200 should be higher than those for haloes that have 
the same numerical value for M v i r . Figure A I shows that this 
is indeed the case. As a result, we generically expect M200 
haloes at the present time to have less substructure than 
haloes with the same numerical value of M V1I (for compari- 
son see Figure 4). 



APPENDIX B: RESCALING MASS 
DEFINITIONS 

In the spherical collapse model, the density within the virial 
radius i? v ir is assumed to be A v i T (z) times the critical den- 
sity at the time of virialization, where A v i r is given by equa- 
tion (1). In principle, this can be used as a criterion to iden- 
tify bound objects. However, because the details involved 
in approximating the timescale to virialization are not well 
understood, this density is not the only one which has been 
used to identify haloes in simulations. For example Diemand 
et al. (2007) and related papers, and Dolag et al. (2004), 
define haloes as being 200 times the mean background den- 
sity. The corresponding radius, -R2uo,f>, is larger than _R v i r at 
z = 0. On the other hand, Springcl et al. (2001); Gao et 



1 1 1 1 1 r 
* 



0.5 - 



O M, 



12 



1 3 



14 



1 5 



log(M) 



Figure A4. Dependence of formation time on definition of halo 
mass: haloes defined to be 200 p C rit(^) will have formed at higher 
redshifts with respect to objects which are A v j r (z) < 200 times 
the critical density. Note, here M corresponds to Af200 i n the case 
of the solid triangles and to Af v i r in the case of the open circles. 



al. (2004); Springel et al. (2005) and related papers define 
haloes as having average densities 200 times the critical den- 
sity. The corresponding radius i?200,c is smaller than _R v i r at 
z = 0. Note that all three definitions become comparable at 
high redshift. 

For each definition, the enclosed mass can be read as: 



Mi 



(Bl) 



where i = {(vir), (200,6), (200, c)}. 

To rescale the different mass definitions we estimate 
R200,b/Rvir and R,200,c/Rvir at various z by computing i?200,c 
and i?200,b from the density profile of our selected haloes 
(which were defined to have sizes R V1I ). We consider 8 simu- 
lation outputs between z = and z — 3. Figure Bl shows the 
evolution of the ratio R200,b/Rvii- The filled pentagons show 
the mean value of this ratio when averaged over haloes of 
all masses; the other data points are for different mass bins, 
here expressed as v = 5\\(z)/S(M). The error bars show 
the rms scatter of the distribution at each z for each halo 
sample. The solid curve corresponds to: 



R 



200,6 



Rv 



= 1.236 



A vir (z) 
A vir (0) 



(B2) 



where A v i r (z) is given by equation (1) in the main text. The 
rms scatter between this curve and the measured points is 

"rmi = 1.1%. 

Figure shows the evolution of the ratio R2oo,c/ Rvii- 
Data points and error bars are like in Figure Bl. The solid 
curve corresponds to: 



-R200,c 
-Rvir 



= 0.746 
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(B3) 
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(1997); this form has a characteristic scale r s , where the 
density profile changes slope, and the concentration is de- 
fined as the ratio of this scale to the virial radius. Thus, 
depending on the density used to define the halo, one might 
have C200 = R200,c/r a (Navarro et al. 1997; Neto et al. 2007; 
Gao et al. 2008), c 2 oo,b = R200,b/r s (Dolag et al. 2004) or 
c v h- = R v [ r /r s (Bullock et al. 2001). Because r 3 itself is a 
physical quantity that does not depend on the halo mass 
definition, the expressions above can be used to transform 
between the different definitions of concentration. 



Figure Bl. Evolution of the ratio i?200,t/^vir- Symbols show 
the simulation data for different masses (here expressed as u); 
solid curve shows equation (B2). 
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Figure B2. Evolution of the ratio i?200,c/ Rvlr- Symbols are as 
in Figure Bl; solid curve shows equation (B3). 



that well fit the mean value computed considering all haloes, 
with a o rms — 0.08%. 

For consistency we check that the ratio i?200,i/^200,c 
approach to unity at high redshift where the universe is mat- 
ter dominated and the dark matter density approaches the 
critical one. 

Equations (B2) and (B3) are also useful for rescaling 
different definitions of halo concentration. Halo density pro- 
files are will fit by the functional form of Navarro et al. 



